package Reports::PsiteReport;

use strict;
use warnings;

my $VERSION = '1.01';
our @EXPORT = qw(gen_psite_report);
use base qw(Exporter);

sub gen_psite_report{
	my %data = %{$_[0]};
	my $report_file = $_[1];
	my $ascore_limit = $_[2];
	my $is_tmt_or_itraq = $_[3];

	my %psites;

	#############	Organizacion de la informacion segun p-sites	################
	#############	          (Llenado del hash psites)	            ################
	foreach my $spectrum(keys %data){
		#Solo reporto los que tengan consenso en el report Global!!!!
		unless(defined $data{$spectrum}{'consensus_mod'}){
			next;
		}
		my @ascores;

		#Solo quiero fosforilados
		my $fosfopeptido;
		if(defined $data{$spectrum}{'ascore_data'}{'peptide'}){
			$fosfopeptido = $data{$spectrum}{'ascore_data'}{'peptide'};
			if(defined $data{$spectrum}{'ascore_data'}{'ascore'}){
				@ascores = @{$data{$spectrum}{'ascore_data'}{'ascore'}};
			}
		}
		else{
			next;
		}

		#Agrupo las isoformas (el array @proteins solo tiene 1 elemento)
		my @proteins;	#=split /\|/,$data{$spectrum}{'consensus_proteins'};
		$proteins[0] = $data{$spectrum}{'consensus_proteins'};

		#El modo en que aqui se agrupan las proteinas NO es optimo, pero por ahora, funciona.
		#La clave de agrupacion es el subconjunto de proteinas asignadas a cada peptido.
		#El problema vendria con dos subconjuntos en que tengamos una interseccion.
		#En ese caso, tendriamos un p-site duplicado para los peptidos que se incluyan
		#dentro de la interseccion.
		#En el dataset "Nov08" hay al menos un caso.
		#O00499	-> O00499|O00499-2|O00499-3|O00499-4|O00499-5|O00499-6|O00499-7|O00499-8|O00499-9|O00499-10|O00499-11->s_282
			#G(214)NK(214)SPS(21)PPDGSPAATPEIR
		#O00499 ->	O00499|O00499-2|O00499-6|O00499-11	-> s_315
			#V(214)NHEPEPAGGATPGATLPK(214)S(21)PSQLR
		#Los valores de ejemplo corresponden a la:
		#Proteina representante->todas las proteinas->aa_#posicion_en_prot_representante
		# y en la segunda linea el peptido correspondiente
		#Es importante tener en cuenta que la localizacion del p-site se da solo
		#para la modificacion dentro de la proteina representante; para las demas
		#"isoformas" deberia obtenerse la localizacion de cada p-site de forma independiente


		my $aa_num;

		#Si agrupo segun proteinas representantes ("isoformas"), este array solo tendra 1 elemento !!!
		foreach(@proteins){
			$_ =~ s/\((\w+)\)//g;
			my ($prot, $aa_num) = ($_, $1);
			unless($aa_num){
				next;
			}
			my @phosphosites = @{phosphosites($fosfopeptido, $aa_num)};

			for(my $i=0; $i<$#phosphosites+1; $i++){
				if(defined $psites{$prot}{$phosphosites[$i]}[0]){
					my %dat;
					$dat{'spectrum'} = $spectrum;
					$dat{'phosphopeptide'} = $fosfopeptido;
					$dat{'valid'} = 0;
					if(@ascores){
						$dat{'ascore'} = $ascores[$i];
					}
					if ($is_tmt_or_itraq eq 'itraq'){
						if(defined $data{$spectrum}{'itraq_info'}){
							$dat{'itraq_info'} = $data{$spectrum}{'itraq_info'};
							my @quant = split /\|/, $data{$spectrum}{'itraq_info'};
							my $ratio;
							if(($quant[0] + $quant[1]) == 0){
								$ratio = 0;
							}
							else{
								$ratio = ($quant[2] + $quant[3]) / ($quant[0] + $quant[1]);
							}
							$dat{'ratio'} = sprintf("%.3f", $ratio);
						}
					}
					elsif ($is_tmt_or_itraq eq 'tmt'){
						if(defined $data{$spectrum}{'tmt_info'}){
							$dat{'tmt_info'} = $data{$spectrum}{'tmt_info'};
							my @quant = split /\|/, $data{$spectrum}{'tmt_info'};
							my $ratio;
							if(($quant[0] + $quant[1]) == 0){
								$ratio = 0;
							}
							else{
								$ratio = ($quant[2] + $quant[3]) / ($quant[0] + $quant[1]);
							}
							$dat{'ratio'} = sprintf("%.3f", $ratio);
						}
					}
					if(defined $data{$spectrum}{'ms2_spectral_info'}{'ms2_file'}){
						$dat{'ms2_info'} = $data{$spectrum}{'ms2_spectral_info'}{'ms2_file'};
					}
					$dat{'proteins_count'} = 1 + $#proteins;
					push @{$psites{$prot}{$phosphosites[$i]}[0]}, \%dat;
				}
				else{
					$psites{$prot}{$phosphosites[$i]}[0][0]{'spectrum'} = $spectrum;
					$psites{$prot}{$phosphosites[$i]}[0][0]{'phosphopeptide'} = $fosfopeptido;
					$psites{$prot}{$phosphosites[$i]}[0][0]{'valid'} = 0;
					if(@ascores){
						$psites{$prot}{$phosphosites[$i]}[0][0]{'ascore'} = $ascores[$i];
					}
					if(defined $data{$spectrum}{'itraq_info'}){
						$psites{$prot}{$phosphosites[$i]}[0][0]{'itraq_info'} = $data{$spectrum}{'itraq_info'};
						my @itraqs = split /\|/, $data{$spectrum}{'itraq_info'};
						my $ratio;
						if(($itraqs[0] + $itraqs[1]) == 0){
							$ratio = 0;
						}
						else{
							$ratio = ($itraqs[2] + $itraqs[3]) / ($itraqs[0] + $itraqs[1]);
						}
						$psites{$prot}{$phosphosites[$i]}[0][0]{'ratio'} = sprintf("%.3f", $ratio);
					}
					if(defined $data{$spectrum}{'ms2_spectral_info'}{'ms2_file'}){
						$psites{$prot}{$phosphosites[$i]}[0][0]{'ms2_info'} = $data{$spectrum}{'ms2_spectral_info'}{'ms2_file'};
					}
					$psites{$prot}{$phosphosites[$i]}[0][0]{'proteins_count'} = 1 + $#proteins;
				}
			}
		}
	}

###################		Para cada p-site, obtengo si es valido para calculos globales	  ##############
########   	(Si sipuera el ascore limite y si no es un espectro MS3 redundante con MS2 valido)##########

	foreach my $prot(keys %psites){
		LAB0:foreach my $psite(keys %{$psites{$prot}}){
			my (@i114, @i115, @i116, @i117);
			my @ratios;
			my $espectros_itraq_usados = 0;
			my @datos = @{$psites{$prot}{$psite}};
			### Si el Ascore es menor a $ascore_limit, no usamos datos itraq para promediar
			LAB1:for(my $i=0; $i < $#{$datos[0]}+1; $i++){
				if(defined $datos[0][$i]{'ascore'}){
					if($datos[0][$i]{'ascore'} < $ascore_limit){
						next LAB1;
					}
				}
				### 1/-Si es un MS3, 2/-y el MS2 esta en el listado, 3/-y el MS2 sale
				###    con ascore>=$ascore_limit, no usamos este MS3
				if(defined $datos[0][$i]{'ms2_info'}){
					for (my $j=0; $j< $#{$datos[0]}+1; $j++){
						if($datos[0][$j]{'spectrum'} eq $datos[0][$i]{'ms2_info'}){  #punto2
							if($datos[0][$j]{'ascore'} >= $ascore_limit){            #punto 3
								next LAB1;
							}
						}
					}
				}
				#Si he llegado hasta aqui, etiqueto el espectro como valido
				$psites{$prot}{$psite}[0][$i]{'valid'} = 1;
			}
		}
	}

####################################################################################################
#####################################################################################################
#print SAL Dumper %psites; #vuelco el hash %psites para debug y control

##########################		Genero el report para los psites	  			#####################
###########   					(calculo medias "on the fly")					#####################

	$report_file = $report_file . '.txt';
	open SAL2, '>', "$report_file";

	foreach my $prot(sort(keys %psites)){
		my $prot_represent = (split /\|/, $prot)[0];
		print SAL2 $prot_represent, "\t", $prot, "\n";
		foreach my $psite(keys %{$psites{$prot}}){
			print SAL2 "\t", $prot, "\t";
			print SAL2 $psite, "\t";
			my @datos = @{$psites{$prot}{$psite}};
			#informacion de promedios (solo con datos validos)
			my (@i114, @i115, @i116, @i117, @ratios_usados);
			my @itraqs;

			for(0..$#{$datos[0]}){
				if($datos[0][$_]{'valid'} eq 1){
					if(defined $datos[0][$_]{'itraq_info'}){
						push @ratios_usados, $datos[0][$_]{'ratio'};
					}
					else{
						#print "espectro ",$datos[0][$i]{'spectrum'}," sin iTRAQ info\n";;
					}
				}
			}

			my $promedio_ratios_verdadero = Math::NumberCruncher::Mean(\@ratios_usados);
			if(defined $promedio_ratios_verdadero){
				$promedio_ratios_verdadero = sprintf("%.2f", $promedio_ratios_verdadero);
			}
			else{
				$promedio_ratios_verdadero = 'NC';
			}

			my $desvest_verdadero = standard_deviation(@ratios_usados);
			unless(defined $desvest_verdadero){
				$desvest_verdadero = 'NC';
			}
			my $CV_ratios_verdadero = '';

			if($promedio_ratios_verdadero eq 'NC'){
				$CV_ratios_verdadero = 'NC';
			}
			elsif($promedio_ratios_verdadero == 0){
				$CV_ratios_verdadero = 'NC';
			}
			elsif($desvest_verdadero eq 'NC'){
				$CV_ratios_verdadero = 'NC';
			}
			else{
				$CV_ratios_verdadero = 100 * $desvest_verdadero / $promedio_ratios_verdadero;
			}

			my $spectra_used = scalar @ratios_usados;

			if(defined $promedio_ratios_verdadero){
				print SAL2 "\t\t\t\t", $promedio_ratios_verdadero, "\t";
				print SAL2 $desvest_verdadero, "\t";
				print SAL2 $CV_ratios_verdadero, "\t";
				print SAL2 "used specs:", $spectra_used;
			}
			else{
				print SAL2 "\t-\t-\t-\t-\t-\t-\t", "used specs:0";
			}

			print SAL2 "\n";

			for(my $i=0; $i < $#{$datos[0]}+1; $i++){
				print SAL2 "\t\t",$datos[0][$i]{'phosphopeptide'}, "\t";
				if(defined $datos[0][$i]{'itraq_info'}){
					my @itraqs = split /\|/, $datos[0][$i]{'itraq_info'};
					foreach(@itraqs){
						print SAL2 $_, "\t";
					}
					if($datos[0][$i]{'valid'} == 1){
						print SAL2 $datos[0][$i]{'ratio'}, "\t";
					}
					else{
						print SAL2 "[", $datos[0][$i]{'ratio'}, "]\t";
					}
				}
				else{
					print SAL2 "\t\t\t\t\t";
				}
				print SAL2 $datos[0][$i]{'spectrum'}, "\t";
				if(defined $datos[0][$i]{'ascore'} ){
					print SAL2 $datos[0][$i]{'ascore'}, "\t";
				}
				else{
					print SAL2 "\t";
				}
				print SAL2 "\n";
			}
			print SAL2 "\n";
		}
	}
	close SAL2;
}


sub phosphosites{
	#unifica fosforilacion y deshidratacion y devuelve una cadena del tipo t_12,
	#donde 12 es la posicion de la modificacion dentro de la proteina
	my ($pep, $aa_num) = @_;
	my @phosphosites;
	$pep =~ s/(\(214\)|\(737\)|\(35\))//g;   #incluyo TMT!
	$pep =~ s/(214\,|\,214)//g;
	$pep =~ s/(737\,|\,737)//g;
	$pep =~ s/(S\(21\)|S\(23\))/s/g;
	$pep =~ s/(T\(21\)|S\(23\))/t/g;
	$pep =~ s/Y\(21\)/y/g;
	$pep =~ s/\(35\)//g;

	my @aas = split //, $pep;
	for(0..$#aas){
		if(($aas[$_] eq 's') or ($aas[$_] eq 't') or ($aas[$_] eq 'y')){
			push @phosphosites, $aas[$_] . '_' . ($_ + $aa_num);
		}
	}
	return \@phosphosites;
}


sub standard_deviation {
	my(@numbers) = @_;
	#Prevent division by 0 error in case you get junk data
	my $divisor = scalar @numbers;
	$divisor = $divisor - 1;
	return undef unless($divisor);
	return undef unless(scalar(@numbers));

	# Step 1, find the mean of the numbers
	my $total1 = 0;
	foreach my $num (@numbers) {
	$total1 += $num;
	}
	my $mean1 = $total1 / (scalar @numbers);

	# Step 2, find the mean of the squares of the differences
	# between each number and the mean
	my $total2 = 0;
	foreach my $num (@numbers) {
	$total2 += ($mean1-$num)**2;
	}

	my $mean2 = $total2 / $divisor;

	# Step 3, standard deviation is the square root of the above mean
	my $std_dev = sqrt($mean2);
	return $std_dev;
}


return 1;          #344

__END__

=head1 NAME

=head1 SYNOPSIS

=head1 DESCRIPTION

=head1 LIMITATIONS

=head1 SEE ALSO

=cut
